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A discrete time model that is capable of replicating the basic features of cardiac cell action 
potentials is suggested. The paper shows how the map-based approaches can be used to design 
highly efficient computational models (algorithms) that enable large-scale simulations and analysis 
of discrete network models of cardiac activity. 
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I. INTRODUCTION 

The heart cells that are directly involved in the dynam- 
ics of its electrical activity include pacemaker and non- 
pacemaker cells. The peacemaker cells generate sponta- 
neous action potentials and are characterized by a slow 
rate of the depolarization. These cells are found in sinoa- 
trial and atrioventricular nodes of the heart and initi- 
ate the propagation of electrical activity throughout the 
heart. The non-pacemaker cells, involved in the propa- 
gation of electrical activity, such as atrial myocytes, ven- 
tricular myocytes and Purkinje cells, have a specific form 
of action potential (AP) which is characterized by five 
main phases including a very rapid depolarization and a 
prolonged plateau pj, see Fig[T] This paper proposes a 
simple computationally efficient model that can replicate 
these phases of AP. 

The non-pacemaker cell has a very negative resting 
potential (phase 4) characterized by the open potassium 
channels and, therefore, K"*" current and the closed fast 
sodium Na"*" and slow calcium channels Ca^"*" . Being de- 
polarized to the threshold voltage of about -70mV the cell 
shows a very rapid depolarization caused by opening fast 
sodium channels Na"*" (phase 0). This phase quickly in- 
creases the membrane potential to positive values where 
the Na"*" channels become inactivated and the dynam- 
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ics of the membrane changes to an initial repolarization 
(phase 1) induced by a short-term transient outward 
current. The plateau of cardiac action potential (phase 2) 
is the result of balanced dynamics between inward Ca^"*" 
current and the outward K"*" current coming through the 
slow delayed rectifier potassium channels. A number of 
other ionic currents are also involved in this phase. As 
the Ca^+ channels start to close while the rectifier K''" 
channels remain open, the membrane potential drops to 
the levels of resting potential forming phase 3. Even from 
this, an overly simplified scenario involved in the forma- 
tion cardiac action potential, it becomes clear that any 
attempt to produce a conductance based model of a non- 
pacemaker cell will lead to a large system of differential 
equations. A number of such models have been proposed 
and are used as accurate models of cardiac cells, see for 
example 0, i, [i| 

The simplified approaches to the modeling of cardiac 
ventricular action potentials include the sets of ODE 
models where each ODE equation phenomenologically 
represents the dynamics of multiple channels (for exam- 
ple, van der Pol equation or the FitzHugh-Negumo 
model [1, 0)- However, due to the high depolarizing 
rates of AP during phase numerical simulations with 
the ODE models require a significant reduction of the 
integration step size which complicates the use of these 
simple models in studies of large scale networks. Here we 
suggest a model of the cardiac AP which is built using a 
discrete time map and designed to significantly increase 
the time step of simulations. A similar approach was 
successfully used before in the design of computationally 
efficient models of spiking-bursting [8,], regular spiking 
and fast and spiking neurons 0, and other neurons. 



II. MAP-BASED MODEL OF CARDIAC AP 

The simplest form of the suggested model is a two- 
dimensional map, which can be written as 



FIG. 1: A sketch of typical action potential of a ventricular 
cell illustrates the main phases associated with different states 
of channel activity. 
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FIG. 2: (Color online) The shape of f(xn,u) and a limit cycle 
generated by subsystem ifTa)) with e{yn) ~ and a fixed value 
of u=-2.55 

where the dynamics of Xn represents the fast changes 
related the phase and the dynamics of j/„ forms the ac- 
tion potential during the remaining phases. The voltage 
of action potential Vn will be defined at the end of this 
section as a liner combination of x„ and y„. 

The function PixmUn) in the right hand side of the 
first equation (fTa|) can be written as 

P{Xn,yn) = (1 - e{yn))f{Xn:Xn-l,u) + eiy„)Xp, (2) 

where the nonlinear function is of the 

form [| 

f{x„,Xn-l,u)= (3) 

a/{l — Xn) + u, if a;„ < 

a + u, if < Xn < a + u and Xn-i < 

— 1, a Xn > a + u or Xn-i > 0. 

Note, that for a description of autonomous dynamics of 
the model the conditions related to the values of Xn~i 
can be omitted. The third argument u will represent a 
linear combination of the function parameter and the 
input variable (e.g. injected current) 

U = f3.^+In. (4) 

Parameter a is a control parameter of the map. The 
dependence of /(x, u) on x computed for fixed values of 
u and e{yn) = 0, i.e. when P{xn,yn) = f{x,u), is shown 
in Figl21 This figure also illustrates a trajectory of an 
uncoupled one-dimensional map ()lap . The limit cycle of 
the map was used in ^] to replicate a sequence of short 
neuronal action potentials - a spike train. Note, that the 
third condition of f{x,u) corresponds to the moment of 
time when Xn reaches its maximum value, i.e. the tip of 
a spike. 

To replicate the action potential of a cardiac ventricu- 
lar cell we shift f{x, u) down where the nonlinear segment 
intersects the diagonal giving birth to stable and unsta- 
ble fixed points. This shift destroys the limit cycle shown 




FIG. 3: (Color online) The dependance of Q{x„, y„) on y„ and 
the trajectory of subsystem (|lb|) started at j/s = 1.3 plotted for 
the parameter values fj, — 0.02 and g — 1.0. 

in Figl2l In addition to that we modulate the 1-D map 
with e{yn), where < e{yn) < 1, see ([2]). One can see 
from ^ that when e{yn) approaches one the x depen- 
dance of function P{x„,yn) flattens and approaches the 
values P(xn,yn) Xp for all a;„. Therefore, the increase 
of e(yn) "deactivates" the motions of the 1-D map (fTa|) 
by forming a superstable fixed point Xp. We will utilize 
this effect to replicate the blocking of Na"*" channels to 
terminate phase 0, see Fig[TJ 

When equation (jlap produces a spike the membrane 
potential rises and the other ionic channels get involved 
the formation of the plateau described above as phases 
1,2 and 3. The dynamics of the cell during this part of 
AP is governed by subsystem (jlb[) . where the right hand 
side can be written in the form 

, , ( ys, if Xn > a + u or Xn-i > 0, 
WKXn^yn)- |g(y^)^ otherwise. 

where 

q{yn) = (1 - M)yn - 52/«(i - ynf- 

Note that the first condition of ([5]) is the same as the 
third condition of ([3]) and corresponds to the moment of 
spike in subsystem (|lap . One can see that at the mo- 
ment of spike the trajectory of the subsystem (|lb[) starts 
at the value y„ = ys and then follows the dynamics of 
the one-dimensional map which is modeled here with a 
polynomial function q{yn)- Function q{yn) is designed 
to form a stable fixed point at yn—0 and a very narrow 
gap between the function and the diagonal where the tra- 
jectories of subsystem (llbp slow down, see Fig. [31 This 
slow motion is used to form the plateau in the shape of 
cardiac action potential, i.e. phase 2. The size of the 
gap and, therefore, the duration of the plateau, is set by 
selecting parameter /z. The other control parameter g of 
the function is used to shape the action potential at the 
transient from plateau to the resting state, i.e. phase 3. 
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The idea behind the selection of these parameter values 
is illustrated in Fig 01 The duration and shape of the 
action potential can also be controlled by the location of 
starting point in equation (|lb[) . i.e. parameter j/j. 
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FIG. 4: (Color online) Waveforms generated by subsystem 
(|lb|) with j/s = 1.3 and different values of parameters fi (top 
panel) and g (bottom panel). 

In order to model the absolute refractory period (ARP) 
that prohibits triggering a new action potential during 
phases 1, 2 and the beginning part of phase 3 we deac- 
tivate subsystem (fTa|) using function e(?/„) which can be 
defined as a step function 



e(2/«) 



if Vn 
if Vn 



< Vth 

> yth- 



(6) 



Here, the value of yth sets the threshold of deactivation. 

Equation ^ closes the feedback between subsystems 
(fla| and (jlbp , and completes the basic part of the map- 
based model. The membrane potential replicated by this 
model can be defined as 

V[mV] = -25mV 4- a;„ x 80mV + y„ x 85mV. (7) 

Examples of APs produced by this model are presented 
in Figini The parameters a and Px of function ([3]), (g]) 
are set to provide a stable fixed point in the subsystem 
(fla| at the negative values of a;„. This takes place when 
the nonlinear part of f{x„) crosses the diagonal. The 
AP was triggered by a pulse of external current /„, see 
(|^. This pulse moves function f{xn,u) up and, if the 



amplitude of the pulse is sufficient, the stable fixed point 
disappears via a saddle-node bifurcation in ^Ta^ . After 
that the trajectory Xn goes up forming a short spike in 
the waveform of x„. The time interval between the trig- 
ger pulse and the spike depends on the amplitude and 
duration of the pulse, see Fig. \5\ The spike initiates 
the motion in subsystem (jlbp by changing its state to 
Un = Us, see ([5]). After y„ has occurred at the high level 
it deactivates subsystem (fTa| by setting e{yn) — 1 and 
keeps it deactivated until ?/„ gets to the levels below yth, 
see dH). 
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FIG. 5: (Color online) Example of ventricular action poten- 
tials computed with formula ([7]). The parameter values are 
a = 3.2, = -2.5780, Xp = -0.8, ys = 1.3, y. = 0.002, 
g = 0.2, yth ~ 0.5 and three amplitudes Ap of triggering 
pulses of external current /„ which is applied at n = 20 for 
one iteration. 



Despite the simplicity of dynamical mechanisms in- 
cluded in the map-based model constructed above the 
model is capable of replicating a number of interesting 
properties of the behavior of ventricular cells. For exam- 
ple it can capture the effects of the coexistence of differ- 
ent oscillation regimes triggered by a periodic sequence 
of pulses. With proper selection of frequency and am- 
plitude of the pulses the model can produce APs with 
the frequency ratio 1:2 or 1:1. These regimes are shown 
in Fig. [S] When in addition to the periodic sequence of 
pulses we add one more pulse unrelated to the periodic 
sequence the result depends on the phase of AP where 
the new pulse occurs. If the additional pulse falls within 
the absolute refractory period it does not affect the oscil- 
lations, see Fig. [5^. If the pulse occurs before or after the 
ARP it can switch the oscillations to the regime with a 
different frequency locking ratio, as it is shown in Fig.[6)D. 
The coexistence of these regimes indicate the presence of 
a memory effect in the dynamics of cardiac AP. In this 
basic model the memory is the result of a transient in 
subsystem (|lal) from the state Xp to the fixed point after 
it has been activated. Similar phenomena caused by a 
memory effect have been studied in [l^, [O, [13 with a 
different type of mapping model. 
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FIG. 6: (Color online) Action potentials triggered by a peri- 
odic sequence of pulses (red circles) and one additional pulse 
indicated by the arrow. The amplitude of pulses A = 0.05. 
The other parameters the same in Fig. [5] (a) Regime of os- 
cillations with frequency ratio 1:2. (b) Transition from 1:2 to 
1:1 frequency ratio caused by the additional pulse. 
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FIG. 7: (Color online) Waveforms of i/„ and AP illustrating 
the effect of electrical restitution modeled by means of ys 
modulation using Eq. [8]with parameter values Sa ~ 0.6, Ss ~ 
0.599, Qs ~ 0.02. The parameter values of the map-models 
are the same as in Fig. [5] except for A — 0.69, g = 0.1 and 
Ai = 0.001. 



III. MODELING OF MEMORY EFFECTS 

Regulation and rate dependance of action potential du- 
ration (APD) is an important property of the ventricular 
cells lillllli]. In the considered map-based model the 
effects related to APD adaptation can be achieved by 
modulation of the subsystem (jlbp by dynamically vary- 
ing the values of parameters /i, g or ys, see ([5]). As an 
example consider the case when the parameter ys in ^ 
is substituted with a new variable = j/s — s„ and s„ 
is described by equation of the form 

_ ( Ss, if Xn > a + u or x„^i > 0, , , 
\'Zs('Sn), otherwise. ^ ' 

where qs(s„) = s„ — gsSn(Sa — Sn)- This equation 
works similar to subsystem (ITb|) . . When subsystem 
(fTa|l generate a spike the trajectory of (|8]) starts at Ss, 
(0 < < Sa) and drifts to s„ = with the rate con- 
trolled by parameter gs, {0 < gs < S^^)- By selecting 
parameter values Sa = 0.6, Ss = 0.599, gs ~ 0.02 one 
can tune the evolution of s„ to replicate the properties 
of the electrical restitution curves. An example of such an 
adaptive behavior is shown in Fig. [71 The figure presents 
a set of waveforms of two consecutive action potentials. 
The first AP is triggered at n=200 while the time of the 
second AP is varied by changing the timing of the second 
trigger pulse. One can see that due to dynamical modu- 
lation of ys the duration of the second AP is significantly 
reduced at the short time intervals between the trigger 
pulses. The APD recovers as the time interval increases. 
The shape of the restitution curve can be controlled by 
the parameters of system ([S]) . 



IV. MODELING OF ELECTRIC COUPLING 

An important component in the simulations of elec- 
trical activity of heart tissue is the model for electrical 
coupling among the cells. This coupling occurs through 
a gap junction between the nearby cells. A number for 
models describing the gap junction dynamics at differ- 
ent levels of so phi stication have been suggested and re- 
ported Here we will consider how the sim- 
plest gap junction model can be used to couple the cells 
in the form of map-based models. We assume that cur- 
rent flowing through the gap junction from cell j into cell 
i is 

I'nji = 9aar>{Vn., - Km) (9) 

where Vn,i is a membrane potential of cell i given by Eq. 
(H) and 9gap is the conductance of the gap junction, i.e. 
coupling strength parameter. 

We will assume that at different phases of AP the cur- 
rent has different effects on the dynamics of the cell. For 
example, at the phases and 4 the effect of coupling 
is more pronounced than during the phases 1,2 and 3. 
Therefore, to insert the gap junction current into the 
map-based model at the phases and 4 we rewrite equa- 
tion dH) by adding /^Jj, i.e. 

where J is the set of nearby cells. In order to capture 
the effect of the gap junction during the phases 1,2 and 
3, when subsystem (jlap is turned off, one can rewrite 
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FIG. 8: (Color online) Waveforms of V„^250 triggered by the 
stimuli (red circles) and images of waves in the circle of 500 
electrically coupled cells. Regime of frequency ratio 1:2 - (a) 
and a change to the regime of 1:1 caused by the proper posi- 
tion of the additional trigger pulse - (b) 



FIG. 9: (Color online) Waveforms of Vn,250 triggered by the 
stimuli (red circles) and images of waves in the circle of 500 
electrically coupled cells, without (a) and with (b) a group of 
"damaged" cells. 



equation (jlbp as follows 



Vn+l 



(11) 



where parameter /Zgap can be selected to set a proper bal- 
ance between the coupling strengths at different phases 
of AP. 

To illustrate some of the effects captured by such a cou- 
pling model we consider a one-dimensional chain with 
periodic boundary conditions which contains 500 cells 
coupled to nearest neighbors. Electrical activity of this 
circle of cells was initiated by excitation of a single cell 
(cell number 250) using a periodic sequence of triggering 
pulses and one additional pulse, whose position in time 
was controlled independently of the periodic sequence. 
The amplitude of the trigger pulses was set A — 0.4 and 
the duration of each pulse was one iteration long. The 
parameters of the cell model were selected as a = 3.2, 



(3x 



-2.5780, 



-0.8, ys = 1.3, /i = 0.001, g = 0.1, 



yth = 0.01, Ss = 0.399, Sa = 0.4 and gs = 0.03. The 
parameters of the gap junction model were set equal to 
ggap = 0.004 and figap = 0.0001. The results of simula- 
tions are shown is Figs. [5] and O 

Figure [5] illustrates the case of the coexistence of two 
regimes of wave activity with the frequency ratios 1:2 
and 1:1. In Fig. [Sja) the additional trigger pulse has a 
phase that inserts perturbations into the wave behavior, 
but docs not change the regime of activity after the tran- 
sient is terminated. When the position of the additional 
pulse was moved closer to the end of the previous AP it 
was sufficient to switch the regime of activity from 1:2 
to 1:1 frequency ratio. This effect is similar to the one 
shown in Fig. [51 but here the memory effects were con- 
trolled directly with ([5]) . The sudden change of frequency 
ratio of the waves is typical for cases of high frequency 
triggering pulses and was demonstrated before in simula- 
tions with more realistic conductance-based models, see 
for example [l^, [l^ . 

Figure [5] shows the effect that can be caused by ex- 
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istence of a small group of "damaged" cells on the wave 
activity of the chain. In this case the period between the 
triggering impulses is larger than the APD of normal cells 
and, when the parameters of all cells are the same, the 
APs with 1:1 frequency ratio is the only stable regime of 
wave activity in the chain. It is illustrated in Fig. [9l^a) 
where the periodic waves perturbed by the additional 
trigger pulse recover very fast. The situation with the 
recovery process changed dramatically after a group of 
30 cells (with indexes i from 150 to 180) was "damaged". 
The " damage" was done by shortening their APD by us- 
ing j/s = 1-2 instead of j/s = 1-3. It is interesting that this 
group of cell did not show signs of "bad' behavior during 
the periodic sequence of triggering pulses, see Fig. [Hb). 
Indeed the three waves at the beginning are almost iden- 
tical to the waves shown in Fig. [D^a) . However after the 
additional trigger pulse perturbed the wave activity of 
the chain the regime of 1:1 oscillations did not recover 
and the system switched to a new high frequency wave 
pattern. This new pattern is characterized by doublets of 
AP waves. The group of "damaged" cells formes a new 
source of wave excitation, see Fig. [SJb) cells from i=150 
to 1=180. 



V. CONCLUSION 

We have considered basic elements of a map-based ap- 
proach to the design of a simple, computationally effi- 
cient model for replication of action potential in a non- 
pacemaker cardiac cell. This paper is focused mainly 
on the basic methods of discrete-time dynamics for the 
model design, rather than on fitting its parameters to 
capture the characteristics of a specific cardiac cell. The 
map-based models tuned for replicating the dynamics of 
specific cardiac cells will be published elsewhere. 
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